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ABSTRACT 


A coarse-grid correction algorithm has been implemented into an implicit upwind Euler 
solver and tested for transonic airfoil problems. The Euler solver uses split-flux formulation and 
pentadiagonal scalar equations, respectively, for the explicit and implicit operators. The multigrid 
sequence starts at the fine grid level, then steps down to each coarse grid level to smooth error 
components using implicit operators. An estimate of residuals can be obtained by two approaches, 
which differ in the level at which the residuals are collected. Both approaches will lead to a work 
reduction factor of 12 for a Mach 0.75 flow at 2° incidence angle on a 65 X 26 grid. The work 
reduction factor is found to increase in proportion to the number of grid levels. 


INTRODUCTION 


The multigrid (MG) method was first advocated in 1980 by Brandt (ref. 1) as a potentially 
efficient technique to solve fluid dynamics equations. Since then, Ni has developed an innovative 
explicit scheme to solve the Euler equations, and Jameson and his coworkers (ref. 3) have attempted 
to combine the MG strategy with explicit as well as with implicit schemes. Jameson's latest work 
has shown an impressive work reduction factor (WRF) as great as 15 for a nonlifting transonic 
airfoil flow on a 128 X 32 grid. These applications share two prominent features; namely, the finite- 
volume formulation and the centered difference approximation to the conventional flux terms. 
Thereby, the error components can be removed and expelled out the computation domain very 
effectively on a succession of grid levels and smoothed by the numerical damping added to the 
difference operators. It is of interest to determine whether the MG strategy will work equally well if 
the flux terms are split into subflux terms according to their eigenvalues and approximated by 
upwind schemes. Because the directional preference of signal propagation might be lost by 
collecting residuals on coarse grids and the numerical damping is generally not required to stabilize 
the results near the shock, the basic mechanisms attributable to the removal of error components 
seem to be missing. To resolve the unsettling issues, the author has incorporated the multigrid 
sequencing into an ADI factorization Euler solver and reported the findings in a 1984 paper (ref 4). 
The MG algorithm is different from other versions in using the corrections at each level rather than 
the residuals at the fine grid level. Although the algorithm is stable, as well as easy to implement, 
and requires no additional storage for intermediate variables, the improvement in WRF is 
relatively small compared to others. It has been found in this study that a much higher WRF can be 
obtained by introducing local time increment and implicit damping to the algorithm. 


SPLIT-FLUX EULER FORMULATION 


The conservation-law form of the inviscid flow equations is given in generalized coordinates £ 

and ip 


« # + F f + G = 0 (1) 

t s n 

Each part of the flux term is composed of three subflux terms which are associated with local 
characteristics A^, Aq. They are of the following form. 
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and y y ij x , q y are treated as local invariants. The three parts corresponding to the characteristic 
values can be cast alternately in the following form. 
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Standard notation is used for flow variables; viz, the density p, the pressure p, and the velocity 
components u and v in Cartesian coordinates, the total internal energy e and (q = u 2 4- u 2 ) /2 , and 
the internal energy e - C V T, which relates to p and p by the equation of state. The sonic speed is 
c = (y p /p) 1/2 ; Y is the ratio of specific heats. The integer f is used to simplify the equation form, in 
which each characteristics component can be identified by t — -1, 0, 1. The conventional matrix of 
characteristics has four components: 


A^ = diag (u — c^, u, u, u + cQ, = diag (U — cq, v, u, v + cq) 

In accordance with the sign of the characteristics, second-order one-sided difference equations are 
used upstream or downstream at each grid point. The order of accuracy is the same on the boundary 
points as those in the field. Details of boundary treatment are given in reference 4 and a 
forthcoming publication. 
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UPWIND FACTORIZATION TECHNIQUE 


Let A v and A u denote the unknown correction vectors abbreviated for (Ap, Au, Ay, Ae) r and 
(Ap, Apu, Apy, Ape) r , respectively. The solution procedure for implicit calculation can be described 
in the following four steps: 


Au. . = (P- I r). . 

UJ \ JUJ 

(/ + - d. 8 U ) A ^ = (r - 1 Ao) w (3) 

I + At8 A -d. 8 )Aw*. . = ls~ l TAw 

on i mj i>j V 

Ay. . = (S Aiy*). . 

uj hj 

where r refers to the residual vector obtained from the explicit operation in equations (1) and (2) 
after replacing the derivatives of flux terms by their difference formulas. The subscripts i and j 
denote the spatial location of a grid point network ranging from / = 1 to imax and j = 2 to j max. 

The operators 8^ and 8^ are for upwind difference and 8^, 8^ are for centered difference associated 
with a damping coefficient, d .. The difference formulas stop switching sides if the magnitude of the 
characteristics is smaller than a prescribed tolerance in order to reduce the truncation errors near 
the stagnation and sonic points. The notation used and the expanded form of the second and third 
equations in equation (3) are given in reference 5. The step increment is generally A£ = CFL • min 
(A£/|A^, Aq/|A n |^) where CFL is the Courant number. It must be equal to or less than 0.25 if 
implicit operations are skipped. A global A t implies that one value of time increment is used for the 
entire domain; a local A t is obtained at each point and used at that point. 



MULTIGRID ALGORITHM 


The implicit procedure for single-level grid calculations may be summarized as follows: 


LAv k+1 = r* . 

IJ lyJ 

£+1 k , A k + 1 

y. . = y. . + Ay. 

1>J IJ lyj 


(4) 


where L is the operator representing equation (3) and r p Au^ j are called the residual and 
correction vectors, respectively. On a coarse grid with spacing twice as large as in all coordinates, 
equation (4) becomes 


LAv 2h =W T r h or 

= < Av h 


(5) 
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On a still coarser grid, 


“"4* = <[(<V r f r ») + %. 

L A V = K | <(< r H - T f H ) + A % 

etc., where W is an operator representing the weighted average of the residuals and corrections on 
the grid points surrounding the coarser grid point, and T is an operator transferring directly the 
values from fine to coarse grid. The correction obtained at each level is interpolated back to the fine 
grid by 


V h = V h + l 2h Ay 2/t 


( 6 ) 


The implicit operator / is a linear function on the boundaries and bilinear inside the domain. The 
weighted operator W serves a very important role to eliminate high frequency modes before passing 
the residuals to coarser grids. Equation (5) is in a general form so that the residuals may be 
obtained from either r or Ao, or from both. For example, letting 


w 2h _ T 2h an( j w 4k _ r 4 h efcc 

The present MG strategy uses the explicit residuals only once on the fine grid level, then uses the 
implicit corrections on other grid levels. A complete cycle of the MG algorithm is illustrated in 
figure 1. 


DISCUSSION 


To help evaluate the performance of the implicit MG algorithm, three parameters (Ap/p )max, 
number of supersonic points NSUP, and stagnation pressure Pt are used for monitoring the 
convergence rate, whereas P t and Mach contours are used for checking solution accuracy. Several 
grids (113 X 34, 97 X 34, 65 X 26, 65 X 22) have beeen considered, but most of the results shown 
here are based on a 65 X 26 grid and for the NACA 0012 airfoil. The baseline calculation uses 
global A t with CFL = 10, as it encounters stability difficulties with local A t even with CFL = 2. 
Figure 2 shows the convergence history and results of a non-lifting case. The magnitude of (Ap/p) in 
the MG calculation is typified by a plateau two orders higher than that in the single-grid (SG) 
calculation. Nevertheless, the MG calculation is at least four times as efficient as the single-grid 
calculation. The WRF is close to 7 if optimal values of CFL = 4 and di = 1 are used instead. Figure 

3 shows the 65 X 26 grid and the Mach contours for a Mach 0.8 flow at 1.25° incidence angle. Figure 

4 is a comparison of the convergence rate between the calculations made from the MG algorithm 
with CFL = 2 local A t and from the single-grid algorithm with CFL = 10. The WRF is estimated to 
be 7 for this case and could have been 13 if CFL = 4 had been used. Note the sharpness of the shock 
above the airfoil and the compression wave underneath. The accuracy of the current results seems 
to be on a par with the conventional SG solution on a 128 X 32 grid. 
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If the MG algorithm uses the explicit residuals on coarse grids; viz, 


L ^2 » = 


the convergence rate is nearly equivalent to that of the current algorithm. Shown in figure 5 is the 
comparison between the two. A comparison of different types of residual collection is shown in 
figure 6. It is seen that the use of local time step increment substantially accelerates the 
convergence rate. 


In the following table the converged values of pressures on the surface and of NSUP are 
listed. The slight differences seen are partly caused by the peculiar behavior of the solution 
technique (eq. (3)), whereby the factorization error affects the converged solutions. 


Parameter 

Pt 

1 Level 

2 Levels 

3 Levels 

Work unit Wu 


1000 

200 

116 

Pmax 

39.46 

40.11 

38.96 

38.72 

Pmin 


15.01 

14.59 

14.78 

NSUP 


86 

88 

88 


M = 0.75, a = 2°; 65 X 26, 33 X 14, 17 X 8 grids 


CONCLUSIONS 


The current MG algorithm is capable of enhancing the stability and the efficiency of the ADI 
factorization technique, and of retaining the excellent shock-capturing capability rendered by the 
split-flux formulation. Satisfactory results are obtained from a 65 X 26 grid for Mach 0.75 flow at 2° 
incidence angle with a work reduction factor equal to 12. The apparent shortcoming of the 
algorithm is in the difficulty of lowering the maximum incremental vector further after a 
reasonable number of iterations. 
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